Adiabatic-Connection-Fluctuation-Dissipation approach to the long-range behavior of 
the exchange-correlation energy at metal surfaces: A numerical study for jellium slabs 
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A still open issue in many-body theory is the asymptotic behavior of the exchange-correlation 
energy and potential in the vacuum region of a metal surface. Here we report a numerical study 
of the position-dependent exchange-correlation energy for jellium slabs, as obtained by combining 
the formally exact adiabatic-connection-fiuctuation-dissipation theorem with either time-dependent 
density-functional theory or an inhomogeneous Singwi-Tosi-Land-Sjolander approach. We find that 
the inclusion of correlation allows to obtain well-converged semi-infinite-jellium results (independent 
of the slab thickness) that exhibit an image-like asymptotic behavior close to the classical image 



potential Vi m (z) 



2 /4 2 . 



PACS numbers: 71.10.Ca, 71.15.Mb, 71.45.Gm 



I. INTRODUCTION 

Two important quantities in the description of a many- 
electron system are (i) the position-dependent exchange- 
correlation (xc) energy per particle e xc (r), which yields 
the total xc energy functional of density-functional the- 
ory (DFT)>i 



E xc [n] = J drn(r)e xc (r), 



(1) 



and (ii) the xc potential V xc (r) entering the Kohn-Sham 
(KS) equation of DFT, which is defined as the functional 
derivative of the xc energy functional: 



V xc (t) 



SE X 



5n(r) 



(2) 



Rigorously, the position-dependent xc energy s xc (r) can 
be obtained as the interaction between an electron at 
r and the coupling-constant averaged charge n xc (r) of 
its xc hole, by using the so-called adiabatic-connection 
formula^— On the other hand, the KS xc potential can 
be obtained by solving the so-called Sham-Schliiter in- 
tegral equation, 5 which relates V xc (r) to the electron 
self-energy of many-body theoryJ^ and which by us- 
ing the Hartree-Fock self-energy reduces to the inte- 
gral equation of the exact-exchange optimized effective 
potential (OEP) schemed Instead, in most of the ex- 
isting electronic-structure calculations these quantities 
have been calculated by invoking the local-density ap- 
proximation (LDA) 8 and its semilocal variants It 
is well known, however, that these approximations fail 
to reproduce the expected image-like asymptotic be- 
haviour of both e xc (r) and V xc (r) at metal surfaces, 11 
an issue of great importance for the interpretation of 
a variety of surface-sensitive experiments such as low- 
energy electron diffraction (LEED)^ scanning tunneling 



microscopy ) ' and inverse and two-photon photoemis- 
sion spectroscopyi 
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The issue of the asymptotic behavior of the xc energy 
and potential at metal surfaces has been the source of 
considerable controversy over the years J£ One of the first 
attempts to describe e xc (z) and V xc (z) at points far out- 
side a metal surface was made by Gunnarsson et al.^ It 
was argued that as the distance z from the metal surface 
grows large enough to make the limiting case of a point 
charge outside a grounded conductor applicable^ (i) 
there are no exchange contributions to order 1/z and (ii) 
correlations are the same as for a classical point charge, 



thus e xc (z -> oo) = V xc (z oo) = V im (z) 



2 /Az. 



The asymptotic behavior of the KS xc potential V xc (z) at 
large distances outside a metal surface was examined by 
Sham (for a semi-infinite jellium) 19 and by Eguiluz at al. 
(for jellium slabs) 20 from the Sham-Shlutcr integral equa- 
tion. Both Shamir and Eguiluz et al2& concluded that 



the Hartree-Fock self-energy yields a 



2 /z 2 behavior 



of the KS exchange potential V x (z) for large z, while the 
correlation self-energy yields the image- like — e 2 /4z limit, 
thereby confirming the long-standing believe that the 
image-potential structure is a pure Coulomb-correlation 
effect. This result disagreed, however, with the semi- 
infinite-jellium calculations reported by Solamatin and 
SahnL— These authors concluded that V x (z) does con- 
tribute to the image-like asymptotic structure of V xc (z); 
they also reached the conclusion that in the asymptotic 
region e x (z) [or, equivalently, half the so-called Slater po- 
tential V x (z)] and V x (z) coincide, but later Nastos^ ar- 
gued that the KS exchange potential should contain the 
entire Slater potential. Most recently, Qian and Sahni 23 
employed the plasmon-pole approximation for the corre- 
lation part of the self-energy to conclude that the asymp- 
tote of the KS xc potential V xc (z) is for metallic densities 
approximately twice as large as the commonly accepted 
— e 2 /4z form and discussed the consequent implications 



2 



of this result on the theory of image states^ 

In an attempt to settle the issue of the long-range 
behavior of exchange and correlation outside a solid, 
fully self-consistent exact- exchange calculations of e x (z) 
and V x (z) have been carried out recently for both jel- 
lium slabs 25 ^ and a semi-infinite jellium i 26 ' 27 For jel- 
lium slabs, it has been proven analytically and numer- 
ically that in the vacuum region far away from the 
surface ef ah (z oo) = V x s (z)/2 = -e 2 /2z^- and 



ySlab^ 



z — > oo) = V x (z) — —e 2 /z 25 which is equiv- 
alent to the well-known results e x (r —> oo) — —e 2 /2r 
and V x (r — > oo) = —e 2 /r that hold in the case of local- 
ized finite systems like atoms and molecules i 17 i 28 For the 
semi-infinite jellium, self-consistent exact-exchange cal- 
culations indicate that (i) the exchange energy per par- 
ticle has an image-like asymptotic behavior of the form 
e x (z —> oo) = —ae 2 /z^ with a density-dependent co- 
efficient a that differs from the jellium-slab a = 1/2 co- 
efficient and coincides with the analytical asymptote ob- 
tained in Ref. I2H but (ii) V x (z) decays as ln(z)/z,— this 
dominant (positive!) contribution not arising from the 
Slater potential V x (z). 

In this paper, we go a step further to report bench- 
mark well-defined self-consistent jellium-slab calculations 
of e xc (z) that include correlation at the level of the 
random-phase approximation (RPA) at least. Our cal- 
culations are based on a combination of the formally- 
exact adiabatic-connection-fluctuation-dissipation theo- 
rem (ACFDT)22 with either time-dependent density- 
functional theory (TDDFT) 30 or an inhomogeneous 
Singwi-Tosi-Land-Sjolander (ISTLS) approach^ As 
pointed out above, the asymptotic forms of the exact- 
exchange energy e x {z) outside jellium slabs and a semi- 
infinite jellium have been found to be qualitatively 
different^ which is due to the fact that for asymptotic 
positions of the electron the exact-exchange (Pauli) hole 
is delocalized and spread throughout the crystal^ 2 - but 
the Coulomb hole screens out the delocalized Pauli hole, 
yielding an xc hole that is localized at the surface! 16 ! 33 
This results in an image-like asymptotic form of the 
position-dependent xc energy e xc (z) of jellium slabs that 
converges quickly to the semi- infinite limit. We find nu- 
merically that this image-like behavior is close to the clas- 
sical image potential Vi m (z) = — e 2 /4z. 



stated otherwise, atomic units are used throughout): 

e xc {z) = i J ' J0j5 J dz ' v ( z ^ z ''^l)n xc (z,z';q), (3) 

where q is a wavevector parallel to the surface, and 
v(z,z';q) and n xc (z, z'\ q) represent the two-dimensional 
(2D) Fourier transforms of the Coulomb interaction 
v(r, r') and the coupling-constant averaged charge 
n xc (r, r') of the xc hole at z' due to the presence of an 
electron at z£~— 

The exact xc-hole charge density, which is related 
to the pair-distribution function and the static struc- 
ture factor of many-body theory, can be obtained from 
the density-response function of linear-response theory 
through the use of the fluctuation-dissipation theorem, 
leading to the so-called adiabatic-connection-fluctuation- 
dissipation formula. For a system that is invariant in two 
directions, we find: 



x dX dujx (z,z';q,iuj) - S[z - z% (4) 
Jo Jo 

where n(z) represents the electron density and 
X X ( Z > z'; q, uj) is the 2D Fourier transform of the interact- 
ing density-response function x A (r,r';u;) at the coupling 
strength A. 

If the interacting X ( z i z '\ Ii entering Eq. (j4|) is 
replaced for all A by the density-response function 
X (z,z';q,oj) of noninteracting electrons moving in the 
effective exact- exchange potential [V e //(z) = Vh(z) + 
V x (z), Vh{z) being the electrostatic Hartree potential] of 
DFT, then Eq. (UJ) reduces to the exact-exchange energy 
per particle: 



£x{z) 



n{z) 



1 -J 



x / dz'ip, l {z,z')(p j (z',z)F lJ {z,z'), (5) 



where ip l (z,z') = £i(z)*£i(z') and 



II. POSITION-DEPENDENT XC ENERGY 

A. Theoretical framework 

Let us consider a jellium slab of width d that is infinite 
in the xy plane (normal to the z axis). The jellium slab is 
invariant under translations in the xy plane, so the xc en- 
ergy per particle at z, defined as the interaction between 
a given electron at z and the coupling-constant averaged 
charge of its xc hole, can be obtained as follows (unless 



dp Jijpk^Jijp&p) 
P vV + {z- z') 2 



(6) 



with J\(x) being the cylindrical Bessel function of first 
order.— Here, £i(z) and &j represent the eigenfunctions 
and eigenvalues of the KS exact-exchange hamiltonian. 
The exact-exchange energy E x is obtained from Eq. ^ 
by replacing e xc (r) with the exact-exchange energy per 
particle e x {z) of Eq. ([5]); and replacing the xc energy E xc 
entering Eq. @ with the exact-exchange energy E x , one 
finds the KS exact-exchange potential V x (z). 
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In the framework of TDDFT, the interacting density- 
response function x A (r, r'; q, uj) obeys the screening inte- 
gral Dyson-like equation^ 

X A (r,r» = x °(r,r',w) + / dndr 2X (r, r', u) 

x{v x {r u r 2 ) + / x A c [n](r 1; r 2 , W )} X A (r 2 , r', u), (7) 

where v (ri,^) = A/|ri — r 2 |, x°( r J r/ 7 w ) is the density- 
response function of noninteracting electrons moving in 
the full KS potential V eff (z) [V eff (z) = V H (z) + V xc (z)\ 
of DFT, and the xc kernel f£ c [n](r,r' ,w) is the func- 
tional derivative of the frequency-dependent xc poten- 
tial of TDDFT. The exact xc kernel is unknown, so 
it has to be approximated. We consider the following 
two approximations: (i) the random-phase approxima- 
tion (RPA), in which the xc kernel (accounting mainly 
for short-range correlation) is simply taken to be zero, 
and (ii) a beyond-RPA approximation, in which the xc 
kernel f^Jn] (z, z'\ q, ui) is borrowed [as in Eq. (43) of 
Ref. l35| from the simple dynamic Constantin-Pitarke 
(CP) uniform-gas xc kernel^ which is accurate in a wide 
range of wave vectors and imaginary frequencies. 

Alternatively, one can follow the recently developed 
ISTLS 31 approach to derive a highly-correlated density- 
response function that has recently proven to yield ac- 
curate xc surface energies^! and an excellent transition 
from three-dimensional to two-dimensional systems*'- 
Within this approach, the density-response function is 
obtained in a self-consistent procedure along with the 
pair-distribution function. 



B. Numerical results 

Figure [T] shows a well-converged RPA calculation of 
Eq. (0} for an electron-density parameter equal to that 
of Al (r s = 2.07) 39 . The corresponding exact-exchange 
semi-infinite-jellium calculation is also plotted for com- 
parison, as obtained from Eq. ([5]) or, equivalently, from 
Eq. (U) by replacing the actual interacting density- 
response function x A ( z : z '; <2S w ) by the corresponding 
noninteracting exchange-only density-response function 
X°(z, z'; q, uS). We note that contrary to the case of exact- 
exchange calculations, where finite-size effects are found 
to be critical^ when RPA correlation is included well- 
converged jellium-slab calculations are a faithful repre- 
sentation of the semi-infinite limit. 

In Fig. we show again the RPA calculation of Fig. [1] 
(thick solid line), but now together with (i) the LDA 
e xc {z) (dashed line) and (ii) the beyond-RPA TDDFT 
calculation that we have performed by constructing a 
jcllium-surface xc kernel from the CP uniform-gas xc 
kernel of Ref. Hll as explained above (thin solid line). 
We have also performed fully self-consistent beyond-RPA 
ISTLS calculations (not plotted in this figure to avoid 
confusion), and we have found a xc-energy curve that 
on the scale of this figure is nearly indistinguishable 
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FIG. 1: (Color online) Exchange-correlation energy per parti- 
cle E X c(z) (in Hartrees) for r s = 2.07, versus z (in units of the 
Fermi wavelength Af). Thick solid line: A faithful represen- 
tation of the RPA xc energy of a semi-infinite metal, as ob- 
tained from jellium-slab calculations. Thin solid line: Exact- 
exchange energy of a semi-infinite metal, as obtained from 
a semi-infinite-jellium calculation. The corresponding fitting 
curves are represented by red dotted lines, as obtained from 
Eq. © with (i) a = 0.30 and z = 0.60 to fit the RPA calcu- 
lation (thick dotted line), and (ii) a — —0.19 and zq — 2.78 
to fit the exact-exchange calculation (thin dotted line). 



from the corresponding beyond-RPA TDDFT calcula- 
tion. Both beyond-RPA (TDDFT and ISTLS) calcula- 
tions are found to capture the correct bulk value (where 
the LDA is exact, by construction, and the RPA is wrong) 
and the correct asymptotics (where the RPA is exact and 
the LDA is wrong). The effect of short-range correlation 
(included in our beyond-RPA calculations) is found to 
be noticeble only in the bulk and in the region near the 
surface, as expected. 

In order to analyze the actual asymptotic (z — > oo) be- 
havior of the (RPA and beyond-RPA) xc energy outside 
the jellium slab, we first write e xc (z) in the image- like 
form 



a 



(z - z )'' 



(8) 



Zq here defining the position of an effective image plane. 
In the case of exact-exchange jellium-slab calculations, 
e x {z) can only be described accurately by Eq. flSJ through 
a fitting of this equation in a vacuum region that is at 
distances from the surface larger than the slab thickness 
d, as discussed in Ref. |26U£ Nevertheless, when correla- 
tion is included Eq. (|8]) nicely reproduces our jellium-slab 
numerical calculations at a few Fermi wavelengths from 
the surface, as shown by the dotted curves of Fig. 1. 

Figs. [3] and 2] show the fitting parameters a and Zq 
that we have found in the RPA for r s — 2.07 as a func- 
tion of the thickness d of the slab. For slab widths that 
are smaller than the distance outside the surface that 
one needs to reach the asymptotic behavior (typically of 
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FIG. 2: (Color online) Exchange-correlation energy per par- 
ticle £ X c{z) (in Hartrees) for r s = 2.07, versus 2 (in units of 
the Fermi wavelength Af). Thick solid line: A faithful rep- 
resentation of the RPA xc energy of a semi-infinite metal, 
as in Fig. [1] Thin solid line: A faithful representation of a 
beyond-RPA TDDFT xc energy of a semi-infinite metal, as 
obtained by constructing a jellium-surface xc kernel from the 
CP uniform-gas xc kernel of Ref. 36. Dashed line: LDA xc en- 
ergy of a semi-infinite jellium. On the scale of this figure, the 
beyond-RPA ISTLS xc energy of a semi-infinite jellium (not 
plotted here) is nearly indisguishable from the corresponding 
beyond-RPA TDDFT calculation (thin solid line). Fitting 
curves are represented by red dotted lines, as obtained from 
Eq. © with (i) a = 0.30 and z = 0.60 to fit the RPA calcu- 
lation (thick dotted line), and (ii) a — 0.30 and zo = 0.30 to 
fit the beyond-RPA TDDFT calculation (thin dotted line). 
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FIG. 3: The parameter a entering Eq. JSJ that fits our jellium- 
slab RPA calculations for r s = 2.07, versus the slab width d 
(in units of the Fermi wavelength Af). 

a few times the Fermi wavelength), we find that both 
a and zq exhibit strong finite-size oscillations. For slab 
widths larger than about three times the Fermi wave- 
length, however, our jellium-slab calculations converge 
nicely to what is expected to be the semi-infinite limit. 
Our converged results for the coefficient a and the 



FIG. 4: The image-plane coordinate zo (in units of the Bohr 
radius) entering Eq. (JSj> that fits our jellium-slab RPA calcu- 
lations for r a — 2.07, versus the slab width d (in units of the 
Fermi wavelength Af). 



TABLE I: The converged a parameter entering Eq. ([8]) that 
fits the asymptotic behavior of the exact-exchange energy and 
the RPA xc energy of a semi- infinite jellium, for various values 
of r s . 



r s 


a x 


n RPA 

XC 


1.5 


0.16 


0.32 


2.07 


0.19 


0.30 


3.0 


0.21 


0.28 


4.0 


0.23 


0.21 


5.0 


0.25 


0.24 


6.0 


0.26 


0.26 



image-plane position zo , as obtained in the RPA and 
beyond RPA (TDDFT and ISTLS), are given in Ta- 
bles I and II for various values of r s at metallic densi- 
ties. The electron-density dependent coefficient a x corre- 
sponding to the semi- infinite exact-exchange e x {z), which 
is included in Table I for comparison, is obtained as 



TABLE II: The converged image-plane coordinate zo (in units 
of the Bohr radius) entering Eq. (JSj that fits the asymptotic 
behavior of the RPA, beyond-RPA TDDFT, and beyond-RPA 
ISTLS xc energies of a semi-infinite jellium, for various values 
of r s . 



r s 


„RPA 
z 


.TDDFT 
z 


z 


1.5 


0.48 


0.23 


0.11 


2.07 


0.60 


0.30 


0.19 


3.0 


0.92 


0.59 


0.38 


4.0 


3.53 


1.16 


1.10 


5.0 


3.02 


1.44 


2.07 


6.0 


2.95 


0.77 


1.29 
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follows^.^ 



k + 2/31n(/3) 
2tt(1 + /3 2 ) : 



(9) 



where f3 stands for the square root of the ratio be- 
tween the Fermi energy ep and the work function W 
(f3 2 = ep/W). There are no entries for beyond- RPA 
a coefficients, since they are very close to their RPA 
counterparts, as expected, which capture the long-range 
behavior of e xc far from the surface into the vacuum. 
We note that the RPA (and beyond-RPA) a coefficient 
agrees (within error bars) for all metallic densities with 
the a = 1/4 classical coefficient for the image potential 
of a classical test charged 

As for the image-plane position exhibited in Table II, 
it is interesting to notice that it is rather sensitive to 
short-range correlations (which are absent in the RPA). 
We also note that the RPA image-plane position zo that 
we have found for r s = 2.07 is close to the image- plane 
position reported in Ref. for V xc (z) and in Ref. for 
the effective surface barrier felt by quasiparticle states 
above the Fermi level of an Al(lll) surface. 

Our calculated image-plane coordinate Zq given in Ta- 
ble II increases smoothly (within the three approxima- 
tions under study) with r s , exhibits a maximum value 
at r s ~ 4, and decreases afterwards. This same pattern 
for zq can be observed for metal spheres from the formula 
/ = W + 1/[2(R + z )]r^ where W is the work function of 
the bulk crystal and R is the radius of the metal sphere. 
Solving this equation for zq and using the values for W, 
R, and I of Table I of Ref. [43], we obtain a pattern for z 
that is similar to that exhibited by our Table [TT1 



III. SUMMARY AND CONCLUSIONS 

In summary, we have reported a numerical study of 
the position-dependent xc energy e xc {z) at metal sur- 
faces, as obtained from jellium slabs, by combining the 
formally exact ACFDT with either TDDFT or an inho- 
mogeneous STLS approach. We have found that the in- 
clusion of correlation allows to obtain (from jellium-slab 
calculations) a faithful representation of the xc energy 
of a semi-infinite system, which exhibits an image-like 
asymptotic behavior of the form of Eq. ([8]) with a coef- 
ficient a that agrees (within error bars) for all metallic 
densities with the a = 1/4 coefficient of the image po- 
tential of a classical test charge. 

The impact of short-range correlation (not present in 
the RPA) has been investigated by either introducing an 
inhomogeneous xc kernel that is borrowed from the sim- 
ple (but accurate) homogeneous CP dynamic xc kernel of 
Ref. |H or following an ISTLS approach. We have found 



that the effect of short-range correlation (included in our 
beyond-RPA calculations) is only noticeable in the region 
near the surface, as expected; this results in a beyond- 
RPA xc energy per particle with an asymptote that also 
has the form of Eq. (|8]) , with an a coefficient that is very 
close to the RPA value (and also to the classical a = 1/4 
value) and an image-plane position zo (see Table II) that 
is sensitive to the introduction of short-range correlation. 

Our self-consistent RPA and beyond-RPA calculations 
lead us to the conclusion that when correlation is in- 
cluded the asymptotics of s xc (z) at metal surfaces agree 
(within error bars) with the classical image potential 
Vi m (z) = —l/Az. The issue of the asymptotic behavior of 
the KS xc potential V xc {z), however, remains unsolved. 
Combining Eqs. flTJ and ©, one easily finds: 



V xc (r) = e xc (r) + / dr'n(r') 



Se xc (r') 
5n(r) 



(10) 



or, alternatively, by using the adiabatic conection for- 
mula to express e^r) as the interaction energy of an 
electron at r with its coupling-constant averaged xc hole, 
one finds: 



V xc {r) = 2e xc (r) + - / drira(ri) ' 



5n(r) 



(11) 

where <?(r, r') represents the coupling-constant averaged 
pair-distribution function. For finite system o ' 28 (and 
also for metal slabs2£), the second term of Eq. (fTTj) [not 
the second term of Eq. (fit)])] does not contribute to the 
asymptotics. However, that is not the case for a semi- 
infinite system, at least when only exchange is taken into 
account. For a semi-infinite metal, the asymptotics of the 
KS exact-exchange potential V x (z) are dominated by the 
second term of Eq. (fTTj). which is always repulsive and 
decays as \n(z)/z£Z- Whether there is a correlation con- 
tribution that asymptotically cancels out this exchange 
term, thus leading to a KS xc potential V xc (z) of the form 
—l/4z [~ e xc (z)] or —l/2z [~ 2e xc (z)], we do not know 
yet. 
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the 7-th and 8-th subbands are partly occupied. We ob- 
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distinguished with the naked eye. 
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